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This paper discusses discrete, linear, time-invariant, nonrecursive, 
finite memory, polynomial smoothing filters for noise that is correlated from 
sample to sample. The wide-sense Markov process is used as a model for 
the noise. Analysis and synthesis of the aforementioned filters are discussed 
in detail and several plots are furnished. A simple method for generating 
discrete, wide-sense Markov noise for simulation is noted. A noise model 
composed of a linear combination of wide-sense Markov processes is de- 
veloped and applied for the case in which the previous model is not suffi- 
ciently accurate. 

I. INTRODUCTION 

A discrete polynomial smoother may be defined in the following 
terms. Consider the random process R(nT), where n is an integer and 
T is the period of the samples at which the process will be of interest.* 
The process will be thought of as comprising a desired component 
R(nT), and a noise component R(nT). It will be assumed that R(nT) 
can be satisfactorily approximated by an rth degree polynomial in nT, 
R{nT). Further, assume that R(nT) is a random process that is wide- 
sense stationary with respect to the sampling instants nT. The fore- 
going situation would occur, for example, in the tracking of a moving 
object whose true position could be represented as an rth degree poly- 
nomial in time, and whose measured position included a random error. 
We will assume that 

E[R(nT)] = (1) 

and denote var [R(nT)] = var [R(nT)] by or, where E is the expected 
value operator of probability theory and "var" indicates "variance of". 

* Symbols used throughout the paper have been collected and defined in a 
glossary (Section IX) for ready reference. 
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Let & B (iT) represent the autocorrelation* function of R(nT), where i 
is an integer. A discrete polynomial smoother of the pth order and rth 
degree is a filter which operates on R(nT) in such a fashion that the 
output C(nT) and the input R{nT) are related by 

E{C{nT)\ = R ip) (nT + r) (2) 

for all w.f Note that the parenthetical superscript (p) denotes "pth 
derivative of the estimate with respect to nT." The quantity r repre- 
sents prediction time; if T is negative, the operation performed is an 
interpolation. 

We will consider linear, time-invariant smoothers which are nonre- 
cursive and have a finite memory. These conditions may be expressed in 
terms of the input-output relationship 

C(nT) = £ W(iT)R[(n - i)T\, ( 3 ) 

t=0 

where the function W(iT) is the weighting function or impulse response. 
Note that W(iT) is defined only at a finite number of points (N points), 
that it is independent of the input (hence the smoother is linear), and 
that it is invariant with the time nT. No previous values of the output 
appear in (3) ; hence the smoother is nonrecursive. The latter restriction 
can often be circumvented, because it is frequently possible to approxi- 
mate a recursive filter by a nonrecursive one. 

The quantity var [C(nT)} = cr c 2 is of interest in two respects. First, 
we may wish to know its value, or better yet, the variance ratio 



or 1 



which is a figure of merit of the smoother. Note that /x 2 is not a function 
of time, since R(nT) was assumed to be wide-sense stationary, and it 
follows that C(nT) is also wide-sense stationary by the time-invariance 
of the smoother. Second, we may wish to find the optimum smoother 
of a class specified by p, r, T, N and T; i.e., we may want to determine 

* In this paper, the term "autocovariance function" will be used to refer to 

E[ZtZ l+T }, 

where Z t represents a zero-mean, wide-sense stationary random process Z evalu- 
ated at time t. "Autocorrelation function" will be used to refer to the normalized 
autocovariance function obtained by dividing the autocovariance function by its 
value at t = 0. 

t The 0th, 1st, and 2nd order smoothers are often referred to as position, ve- 
locity, and acceleration smoothers. 
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the weighting function W(iT) which yields the minimum value of 
(Tc (or ix 2 ) under the preceding conditions. 

If U(nT) has an autocorrelation function of arbitrary form, it may 
be shown, using (1), (3), and (4), that 

,=0 >=0 

where IF, = W(iT) and W j = W(jT). In general, this is a compli- 
cated expression. In previous treatments 2,3,4,5 of discrete polynomial 
smoothers, simplification of (5) has been achieved by assuming that 
the power density spectrum of the noise component of the input is 
white, so that 

fl (» = j) 

*«[(i-J)T]= \ (6) 

This yields the simple form 

M f -i;V\ (7) 

»'=0 

However, the assumption that the noise is uncorrelated from sample 
to sample is not justified for many physical systems because the noise is 
restricted in its rate of change. This is particularly true for mechanical 
and electromechanical systems. It will be shown that correlated noise 
may be represented by the wide-sense Markov process as a first-order 
approximation, or by a linear combination of such processes as a better 
approximation, with appreciable simplification of (5) still being obtained. 
By "represent" we refer to the approximation of one autocorrelation 
function or power density spectrum by another. In discussing smoothers, 
our primary interest is in the behavior of the generalized second moment 
of random processes, and further delineation of the character of these 
processes is not necessary. 

II. WIDE-SENSE MARKOV NOISE MODEL 

A rigorous definition for the wide-sense Markov process may be found 
in Doob. 6 It will be sufficient for our purposes to characterize the wide- 
sense Markov process in an alternative fashion, which Doob 7 has shown 
to be equivalent to the original definition. A wide-sense stationary, con- 
tinuous random process will be called wide-sense Markov if it has the 
autocorrelation function 
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<p(r) = exp(- fir), 



^ 0. 



(8) 



The quantity ft will be called the "noise bandwidth." By using the even- 
ness property for autocorrelation functions of real, wide- sense stationary 
random processes, (8) may be written as 



$(t) = exp(- ft | t|). 



(9) 



If a wide-sense Markov random process is real and Gaussian and has 
zero mean, then it is also strict-sense Markov. The strict-sense Markov 
process is denned as a random process for which 



Pr[F(U ^ X | Y(h), 



F(f„_x)] = Pr[Y(t„) =g X | Yitn-J] (10) 



with probability 1 for each X, all t x < ■■• < t„ , and all n. We may 
say in an intuitive manner that a strict-sense Markov process is a proc- 
ess with a structure such that any value of the process is directly related 
only to the immediately preceding value. 

One might consider higher-order Markov processes ("related" to 
several preceding values) as a better approximation for correlated noise, 
but it appears that using a linear combination of the simple wide-sense 
Markov processes gives a more manageable expression for n . 

For a discrete wide-sense Markov process with equally-spaced samples, 
we may write the autocorrelation function as 



*(t) = exp(-ft | T\)Cb r (T) 
where Cb r is the comb function defined by 

CO 

Cb T (r) = E*(t- iT). 



(ID 



(12) 



S(w) 




Fig. 1 — Baseband component of normalized power density spectrum for dis- 
crete wide-sense Markov process. 
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Fig. 2 — Control system used in evaluation of wide-sense Markov noise model. 

The normalized f power density spectrum, obtained by Fourier trans- 
formation of (11), is 

8(f) = 



2V 



(2tt/) 2 + ft 2 T 



*7nCb llT (f), 



(13) 



where * indicates convolution. The baseband component of this nor- 
malized power density spectrum is illustrated in Fig. 1. Note that the 
half-power point occurs at / = 9,/2-ir. 

Use of the wide-sense Markov noise model reduces (5) to 

M 2 = Y,T,WtW,a x *-\ (14) 



,=o y=o 



where 



a = exp(-flr) 



(15) 



and is called the "intersample correlation." For some weighting func- 
tions, (14) can be simplified much further by evaluating the sums, using 
the finite difference calculus. 

As one illustration of the improvement in accuracy obtained by repre- 
senting correlated noise as wide-sense Markov rather than white, con- 
sider the control system of Fig. 2. White noise is filtered by the continu- 
ous system such that the normalized power density spectrum at the 
input to the sampler becomes 



S(u) = 



G(i«) 



where 



= i£|<?0»)f*o = 4.79. 



(16) 



(17) 



t Normalized in the sense that this is the Fourier transform of the autocorrela- 
tion function. The power density spectrum is usually denned as the Fourier trans- 
form of the autocovariance function. The normalized power density spectrum is 
equal to the power density spectrum divided by the variance. 
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Fig. 3 ■ — Normalized power density spectra. 



Hence 



SM = 



12.6(o 2 + 4) 
a) 4 + 29.3o> 2 + 241.5 



(18) 



We can fit models to the true noise process as if all processes were con- 
tinuous, and following this, introduce the sampling operation. The 
output-input noise variance ratio \? of the digital system has been 
computed for the case of a first-order cascaded simple averages smoother f 
with the following weighting coefficients: 

f 0.028257 (0 ^ i ^ 11) 

Wi = I (12 ^ i ^ 23) (19) 

-0.028257 (24 ^ i ^ 35). 

The true noise process has i* = 0.0376. Use of the wide-sense Markov 
model yields n 2 = 0.0339, while use of the white noise model yields 
ii = 0.0192. 

f See Section V for the definition of this smoother. 
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Fig. 4 — Autocorrelation functions. 

The normalized power density spectra and autocorrelation functions 
of the true noise process and the wide-sense Markov model are illustrated 
in Figs. 3 and 4. The parameter fi has been picked equal to the half- 
power point of the power density spectrum of the true noise process, 
9.68. 



III. MOMENTS OF THE WEIGHTING FUNCTION 

The moments of the weighting function of a smoother are important 
characteristics, since the requirement (2) which specifies the desired 
output of the smoother is conveniently expressed in terms of them. 
The moments will be useful in comparing smoothers for equivalence as 
to meeting (2), and in determining the optimum weighting function for 
a class of smoothers. The </th moment M q of the weighting function 
will be denned as 



M q = E (iT) 9 Wi . 



(20) 



To express (2) in terms of moments, we proceed as follows. Substi- 
tuting (3) and (1) in (2) we obtain 



£ W { R[(n - i)T] = R {p \nT + r). (21) 

i=0 

Now R(t) will be approximated by R(t), which may be expressed in the 
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Taylor series form 



R( t ) = ± fi9( ; ir) (t - nT)\ (22) 



Substituting (22) in both sides of (21) and rearranging, we obtain 

(23) 



t (-Q'iW) g (iTyWi = £ pnT) r ,_ p 
Po ql »=o i=v (q - V) ! 



Considering (23) term by term, and using (20), we obtain 
f0 (0 ^ q < p) 

(-l)V (q = v) 

(~ lY t r 9 " p (p < q ^ r). 

. (q - v) I 



M g = «! 



(24) 



It should be noted that the weighting function obviously has moments 
greater than the rth; however, the condition (2) does not fix their values. 



IV. OPTIMUM SMOOTHERS 

By "optimum smoother" we mean that smoother of the class specified 
by p, r, r, N, and T whose weighting function yields the minimum pos- 
sible value of ix 2 . Optimum smoothers are often not implemented because 
of the amount of storage and computation required. However, they 
provide a standard of comparison for the systems that are implemented. 

To find the weighting function of the optimum smoother of a class, 
the quantity \i is minimized under the constraints (24), using La- 
grange's method of undetermined multipliers. Blackmail 8 has carried 
out the minimization in matrix form for a general input noise process 
(any autocorrelation function). The optimum smoother is specified by 
the matrix equation 

W = P~ l A (1P- 1 A )~ l M , (25) 

where ~ indicates "matrix transpose." The variance ratio n for the 
optimum smoother is given by 

m 2 = M(AP~\A )~ l M. (26) 
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The matrix W is a column matrix representing the weighting function 
at the points t = iT, i.e., 



W = 



W 

Ww-i. 

P is the autocorrelation matrix of the input noise process, 

"1 *r(T) *«(2D 

*a(T) 1 # a (D 



(27) 



P = 



**(2T) 



l 

**(r) 



i 



_**[(# - 1)7*] #«[(JV - 2)rj #«[(iV - S)? 1 ] 
^4 is the "age" matrix, 



A = 



_1 (JV - 1)7' [(iV - l)7f • 
and J\f is the column matrix of moments, 

'Mo 



J/ = 



Mi 



**[(N - 


- DTI"] 


*n[W - 


- 2)7'] 


**[(# - 


-3)7-] 


1 





; (28) 



1 





... o 


1 T 


T 2 


• • • r 


1 27' 


(27') 2 


••• (25T) r 


1 37' 


(3T) 2 


••• (STY 



[(N - 1)7T_ 



; (29) 



(30) 



Unfortunately, (25) and (26) are very difficult to evaluate literally 
except in the simplest cases. However, they can be evaluated numerically 
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by a digital computer. The inverse of the autocorrelation matrix for 
wide-sense Markov noise is readily determined to be, in literal form, 





P-i = 



1 - a 2 











-a 1 + a 2 -a 

-a 1 + a 2 -a. 

-a 1 + a 2 



1 + 



a' — a 



-a 1 



• (31) 



The principal operation, aside from the matrix multiplications, is the 
inversion of the (r + 1) X (r + 1) matrix AP~ l A. 

Blackman 8 has evaluated (25) and (26), assuming that the noise is 
wide-sense Markov, for zero prediction time smoothers with p = 0, 
r = and p = 1, r = 1. For the former, 



1 



Wi = I 



N - (N - 2) a 



1 - 



N - (N - 2) a 



(i = 0, N - 1) 



(i= 1,2, ...,JV-2) 



(32) 



and 



For the latter, 



2 

M = 



1 + a 



Wi = 



N - (N - 2)a' 

(1 + ^[1 + V (N - 2)] 



(33) 



T (N - 1) {[1 + i(tf - 1)][2 + *(tf - 1)] + [1 - >? 2 ]} 

(»"-0) 



! (AT - 1 - 2i) 



6 

r (AT - 1) {[1 + ^(AT - i)][2 + v(N - 1)] + [1 - v 2 )} (34) 

(« - 1,2, ...,Af-2) 

3 (1 + iy)[l + ^(A^ ~ 2)] 

T (N - 1) {[1 + >>(# - 1)112 + v(N - 1)] + [1 - ?)) 

(i = N - 1) 
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and 



M = 



12t7 



T n - (N - 1){[1 + v(N - 1)][2 + n(N - 1)] + [1 - t]} ' 



where 



1 - a 
1 + a 



(35) 



(3G) 



These optimum weighting functions and variance ratios have been 
plotted in a normalized form in Figs. 5, 6, 7, and 8. The ordinates for 
the 1st order, 1st degree smoother are given in terms of the smoothing 
interval T s = (N — 1)T. The curves are plotted for the parameter 
B = SlT s , which may be thought of as a noise-smoother "bandwidth 
ratio." The asymptotes for the above curves, as JV — > <» (with T a 
and ft fixed), are derived in Appendix A. 

Let us consider the behavior of these curves from a physical view- 
point. For wide-sense Markov noise, the noise autocorrelation function 
is positive and monotonically decreasing with time. Hence, if the num- 
ber of samples smoothed, N, is increased with the smoothing interval 
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Fig. 5— Optimum weighting function: Oth order, Oth degree smoother (r = 0, 
n = 6). 
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Fig. 6 — Optimum weighting function : 1st order, 1st degree smoother (r = 0, 
n = 6). 

T s and the noise characteristics remaining - fixed, the inters-ample 
correlation will increase. Although each additional sample provided to 
the smoother gives additional information, the information added 
eventually approaches zero due to the increasing correlation. Now a 
smoother can reduce its variance ratio only by obtaining more informa- 
tion about the noise or by making better use of the information it 
already has. An optimum smoother makes the best use of the informa- 
tion available to it. Consequently, the variance ratio of an optimum 
smoother operating on a signal which includes wide-sense Markov 
noise (or any noise whose autocorrelation function is positive and de- 
creases monotonically with time) must approach a constant as N in- 
creases. 



V. CASCADED SIMPLE AVERAGES SMOOTHERS 

Cascaded simple averages smoothers are a class of smoothers developed 
by R. B. Blackmail. 1 A cascaded simple averages smoother of sth order 



SMOOTHING FILTERS 



2133 



0.22 



0.18 



H-- 



0.08 



0.06 



0.02 



































I 






























V 






























\ 






















PTO 
\ 
\ 


rES 








V 


^__^ 












R- 








































































































































































































































20 












\v 






























$ 
















































30 










































40 


























50 





























fiO 






































7 


3^ 8 


7* 


r~: 





































































1 















4 





C 





e 





■ 






Fig. 7 — Noise variance ratio: optimum Oth order, Otli degree smoother. 

approximates an optimum (with respect to white noise) sth order, sth 
degree, zero prediction time smoother. It may also be used to approxi- 
mate smoothers that have been optimized with respect to wide-sense 
Markov noise. The approximation involves using only the values K, 
—K, and for the weighting coefficients, where K is some constant. This 
smoothing method reduces the amount of storage and the number of 
arithmetic operations required, at the cost of a slight increase in n 2 
over the optimum method. 

The weighting functions of cascaded simple averages smoothers of 
Oth, 1st, and 2nd orders are as follows (respectively) : 



W< = i 

1 N* 



(37) 
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Fig. 8 — Normalized noise variance ratio: optimum 1st order, 1st degree smoother 



Wi = 



N 2 T f 



4MN-1) ( 0SiS f-l) 



(38) 



*& ~ *> ( 2 -N<i<N-l) 

WT a \3 - = V 



and 



TT< = i 
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N 3 T S °- 



36 (JV - 1)'' 

JV 3 2V 



36(iV - if 

N 3 T a * 



(.*<*S-i) 



.6 - - 3 



(|iV ^i^N - lY 



(39) 
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Fig. 9 — Weighting function for 1st order cascaded simple averages smoother. 

where N is a multiple of 3 in (38) and a multiple of 6 in (39) . The weight- 
ing functions for 1st and 2nd order smoothers are plotted in Figs. 9 
and 10, respectively. 

The variance ratios for 0th, 1st, and 2nd order cascaded simple 
averages smoothers for a wide-sense Markov noise input are, respec- 
tively : 



-, JV— 1 AT— 1 

' = \ iM mr - T £ ^ 8 »" W < sgn *'«"""■ 



and 



where 



2 [36(^-1)7^^ w w hW | 

L A' d is- J i-o >-o 

f-i (Wi < o) 

syn Wi = < (TT< = 0) 
1 (Wt > 0) . 



(40) 
(41) 

(42) 



(43) 



By use of the finite difference calculus, (40), (41), and (42) may be 
simplified to 
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Fig. 10 — Weighting function for 2nd order cascaded simple averages smoother. 



M = 



1 + a , 2a(a N - 1, 



+ 



= 1 



(S)'W 
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*\3AT(1 - a) 



/ N n 2N/S AT/3 | n \ 
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[NO - «)f~ 



(44) 



(45) 



and 



= 2 \~ 



s)'( 



A T- 1 



1 4- 



62V/6 



+ 



<(a* - 2a""" - « w/ ' + 2a A/2 + 3a y/ ' - 3) 
[NO ~ «)]* 



(46) 



[3i\T(l - a) 
respectively. 

In Figs. 11, 12, and 13, the variance ratios have been plotted in 
normalized form for 0th, 1st, and 2nd order cascaded simple averages 
weighting functions, respectively. The ordinates are n , T s n , and 
Tap, respectively. The curves are plotted for the noise-smoother 
"bandwidth ratio" B = &T 3 . The asymptotes for the above smoothers 
as N — > oo (with T s and Q fixed) are derived in Appendix A. Note that 
the expressions simplify appreciably for larger values of B, the expo- 
nential terms becoming negligible. 

The behavior of these variance ratio curves is somewhat different 
from those for the optimum smoother. They do not necessarily decrease 
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Fig. 11 — Noise variance ratio: Oth order cascaded simple averages smoother. 

monotonically with N, even though they have asymptotes similar to 
the optimum curves. This is due to the fact that the smoothers are not 
optimum, and therefore the information about the noise is not neces- 
sarily utilized in the best manner. Consequently, as N increases, change 
in variance ratio may be due to changes in the utilization of the infor- 
mation available as well as changes in the information available, and 
the change cannot be readily predicted. 

Note that the curves for all three orders of smoothers (Figs. 11, 12, 
and 13) either have a minimum at some finite value of N or approach 
a minimum as N — > <» . These minima are more or less broad. In specify- 
ing a smoother, it is advantageous to choose the lowest value of N for 
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Fig. 12 — Normalized noise variance ratio: 1st order cascaded simple averages 
smoother. 

which /x 2 is reasonably close to the minimum. Note that the neighborhood 
of the minimum variance ratio as a function of N is reached at lower 
values of N as B decreases (intersample correlation a increases for 
fixed T a ). This is reasonable physically, since the value of smoothing 
a larger number of samples decreases as these samples become more 
highly correlated. 



VI. SYNTHESIS OF POLYNOMIAL SMOOTHERS 

In general, the polynomial smoothers we have been discussing are 
classified by the parameters p, r, V, N, and T.f It would be convenient 



t The optimum smoother is also classified by the parameter a. 
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Fig. 13 — Normalized noise variance ratio: 2nd order cascaded simple averages 
smoother. 

to be able to synthesize the smoother in terms of sth order, sth degree, 
zero prediction time components, where p ^ s ^ r. Note that the com- 
ponents are functions of s, N, and T only; hence their characteristics 
could be specified fairly simply. Further, several smoothers with different 
parameters p, r, and V but the same N and T could be synthesized with 
common components by weighting these components differently. Finally, 
the above breakdown permits any polynomial smoother of the class 
considered in this paper to be constructed from cascaded simple averages 
components. The derivation and procedures discussed in this section are 
valid for discrete polynomial smoothers in general and are not restricted 
to optimum smoothers or to particular input noise power density spec- 
tra. 

Consider the linear combination of sth order, sth degree, zero predic- 
tion time components shown in Fig. 14. Let W»i represent the value of 
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Fig. 14 — Synthesis of pth order, rth degree smoother from sth order, sth degree 
components. 

the weighting function of the .sth component at the sample with age 
IT. Let M aq be the 5th moment of the weighting function of the sth 
component. From Fig. 14 it will be seen that the "over-all" weighting 
function Wi of the entire smoother is related to the component weighting 
functions by 



Wi = £ K a W si . 

■H=p 

Now, using (47) and (20), 

M q = T" £ i q Wi = T" £ K s t, i"W,i = E Wf. 



£ 

itmO 

From (24) we obtain 



j=p t=o 



M sq = 



U-l) q ql (s = q) 
JO (s > q). 

Substituting (49) in (48) we get 



8-1 



M q = 2 KM. g + K q ( -l)*ql 



(47) 



(48) 



(49) 



(50) 



If the linear combination of components is to be equivalent to the pth 
order, rth degree smoother, then (24) must be satisfied. It follows that 
we must have 



K. = 
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(s = p) 



( —^ Z K u M ua (p<s^ r). 



(51) 



(s — p)\ s! 

It should be noted that a linear combination of optimum components 
will not necessarily be optimum unless the outputs of the components 
at a common time are uncorrected. 

The synthesis procedure proper consists of finding a smoother or set 
of component smoothers which produces the desired output (2) with 
the least total error e compatible with a simple implementation. In the 
case of recursive smoothers, stability must be considered; the latter 
topic is adequately covered in standard texts on control theory. The 
total error e is given by 

(52) 



r - i M 
e = \(Tc T «rl] 



where e T is the truncation error 

e T = R ip) (nT + T) - R (p) (nT + r). (53) 

Alternatively, using (21), we may write (53) as 

e T = R (p) (nT + r) - E WiR[(n - i)T\. (54) 

»=o 

Synthesis involves the choice of type of filter (optimum, cascaded 
simple averages, etc.) and the selection of r, N and T. For convenience, 
the parameters will be selected in the alternate form r, N, and T a . 

The selection of r is based on the requirement that r ^ p and the 
direction of change in e as r increases. Now <r c mcreases and t T decreases 
(in general) with increasing r. The rate of increase of <j c with r is such 
that smoothers with r > 2 are seldom used in practice. 

The selection of iV and T„ will be a trial-and-error process based on 
achieving a near minimum in e while keeping N as small as possible 
(for simpler implementation). In the case of a set of components, each 
component may have a different value of T a provided the values of T 
are the same. Figs. 7, 8, 11, 12, and 13 will be useful in calculating <r c 2 . 
When calculating the over-all output noise of a set of component smooth- 
ers, it will be useful to know that the noise outputs of Oth, 1st, and 2nd 
order cascaded simple averages smoothers are all mutually uncorre- 
lated, though this is not true for all orders. 

The problems involved in estimating truncation error have been 
discussed by Hamming 10 in some detail. We will make the simplifying 
assumption that the truncation error e T of an rth degree smoother may 
be approximated by using (54) with R(t) considered as an (r + l)th 
degree polynomial. Thus R(t) may be expressed 
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r+1 d(9> 



R( t ) =Z R9( f } (t-nT)\ (55) 

«=o ql 



Substituting (55) into (54), and using (22) and (26), we obtain 



_(r-p + l)! (r + 

Blackman 11 has calculated the (r + l)th moments of Oth, 1st, and 
2nd order cascaded simple averages smoothers as \T S , —T 8 , and ST S , 
respectively. Hence, the truncation errors for these smoothers may be 
calculated from (56) as $T.R {r+l) (nT). 

VII. GENERATION OF DISCRETE WIDE-SENSE MARKOV NOISE FOR SIMU- 
LATION 

It is frequently desired to simulate the performance of discrete smooth- 
ing filters and perhaps larger discrete systems of which they may be a 
part. Standard techniques are available for simulating discrete white 
noise by generation of a sequence of uncorrected pseudo-random num- 
bers. 1213 It is relatively easy to generate discrete wide-sense Markov 
noise from such a sequence, due to the simple correlation structure of 
the wide-sense Markov process. The foregoing is another advantage in 
using the wide-sense Markov model to represent correlated noise. 

Let { F„} be the desired discrete wide-sense Markov noise and {X n } 
be a sequence of uncorrelated random numbers of zero mean and unit 
variance. Then Y» may be generated as 

F x = aX x , (57) 

F„ = «F„_i + aVl - « 2 X n , (n > 1), (58) 

where a 2 is the variance and a the intereample correlation of the wide- 
sense Markov noise. 

Since F„ is in effect a linear combination of the X n _,- , i = 0, • • • , 
n — 1, it follows that if the X„_,- are jointly Gaussian, then the F„ are 
jointly Gaussian. 

VIII. DISCRETE SMOOTHING FILTERS BASED ON A MODEL USING A LINEAR 
COMBINATION OF WIDE-SENSE MARKOV PROCESSES 

In some cases the simple wide-sense Markov noise model may not be 
a sufficiently accurate representation of a physical noise process. A 
better model may be obtained by approximating the known or assumed 
noise process by a linear combination of wide-sense Markov noise proc- 
esses. We may approximate the autocorrelation function by wide-sense 
Markov autocorrelation functions, or, equivalently, we may approxi- 
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mate the normalized power density spectrum by wide-sense Markov 
normalized power density spectra. For the purpose of making the pre- 
ceding approximations, we can work with the process as though it were 
continuous, later introducing the sampling operation. In a parallel to 
the use of Fourier series to analyze the behavior of a complicated wave- 
form in a linear system, the wide-sense Markov autocorrelation func- 
tions may be used to analyze the behavior of a complicated correlated 
random noise process in a linear discrete system, by applying the prin- 
ciple of superposition. It is possible to synthesize discrete smoothers 
using this more complex model. 

Further, discrete random noise of arbitrary power density spectrum 
may be generated in an approximate manner for simulation purposes by 
a suitable linear combination of wide-sense Markov noise components. In 
the preceding applications, the use of the wide-sense Markov noise com- 
ponents is simpler and more efficient than use of the actual noise process. 

There are two types of approximations that can be made. One is a 
cut-and-try type of approximation in which one tries various linear 
combinations of wide-sense Markov noise components with the band- 
widths of the components not necessarily being integral multiples of 
some fundamental bandwidth. The other approach is to use a linear com- 
bination of orthonormal functions of wide-sense Markov components. In 
the latter approach, the band widths of the components are integral 
multiples of a fundamental component. In either case, we may write 

*(t) = £4,exp(-«,|r|) (59) 

or 

SM=t,A v -^- 2 . (60) 

„=i \l v * + CO 2 

Note that the sum of the coefficients A „ must be equal to 1 . In the ortho- 
normal approximation, 

a, = vu (61) 

and the A v will have a definite form. This is shown in the following sec- 
tion. 

8.1 Orthonormal Approximation 

Laning and Battin 14 and Lee 15 have developed orthonormal approxi- 
mations for an arbitrary autocorrelation function and an arbitrary 
normalized power density spectrum. These approximations are in terms 
of components which will be recognized as wide-sense Markov auto- 
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Table I — Values of Coefficients ct„ 



k 


Chi 


Ckl 


Ck3 


Cki 


CkS 


1 


1 










2 


2 


-3 








3 


3 


-12 


10 






4 


4 


-30 


60 


- 35 




5 


5 


-60 


210 


-280 


126 



correlation functions and normalized power spectra, respectively. We 
shall develop the approximation in somewhat different form. 
The set of functions 



I 



**(»■) = Z) c <» v^fi exp ( —va I t I ) 



(62) 



can be made orthonormal on the interval — °o < T < » by proper 
choice of the coefficients c* P . These coefficients are listed in Table I 
for values of k up to 5. 

These functions may be used to form an orthogonal expansion of any 
piecewise continuous even function (and hence any piecewise continuous 
autocorrelation function) on the interval - oo < T < ». We may write 



where 



*(t) = 53 <W**( r )i 



a*. = / *(t)**(t) dr. 

J—zc 



(63) 



(64) 



If we take z terms of the series expansion and denote the corresponding 
partial sum <S(t), we may group terms to obtain 



z 

<p(r) «$ (r) = 2 Av exp (-1*2 | r | ), 



i'=i 



where 



A v = 2 a * c *» Vfcfl. 



a-=i> 



(65) 



(66) 



The coefficients A„ for 2^5 are given by (note that if z < 5, then 
0* = for k > z) 

A, = Vn k + 2V2 a 2 + 3\/3 a a 4- 4y/l a 4 + 5\/5 aj, (67) 

^ 2 = -3\/ft [a/2 02 4- 4\/3 a 3 + 10\/4 04 4- 20V5 a B ], (68) 
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At = lQ\/a Wl «a + 6\/4 a 4 + 21\/5 05], (69) 

A A = — 35-\/« [V4 fl 4 + 8\/5 oj, (70) 

.4 B = 126\/fl[\/5a B ]. (71) 

Now let S k (co) bo the Fourier transform of $a-(t). 
Then 

&(«) = £ Ct ,Vfcfi , f,"°_ a . (72) 

r=i (viz;- -+■ co- 

It can be shown, using Parseval's theorem, that the set of functions 
{(1/V / 2t)*S'a.(co)J is orthonormal on the interval — oo < co < ». Hence 
we may expand any piecewise continuous even function (and thus any 
piecewise continuous normalized power density spectrum) on this 
interval. We may write 

S(co) = 2 «*&(«), (73) 

where 

a, = -L f S(«)&(«) dco. (74) 

From Parseval's theorem it will be seen that the a k in (74) are the same 
as those in (C>4). If we take z terms of the series expansion and denote 
the corresponding partial sum a(«), we may group terms as before 
to obtain 

8M m 6M = t A. , ffl 5 , (75) 

r=i (yi2J ,! ■+• ur 

where the .4,, are given by (67) through (71). Note that (65) and (75) 
form a Fourier transform pair. Thus, if we have approximated a noise 
process in terms of normalized power density spectra or autocorrelation 
functions, the alternative approximation can be immediately obtained. 
The quantity rft represents the half-power point for each wide-sense 
Markov component. 

Simple rules for the selection of ft for a particular expansion cannot be 
established; it is a matter of judgment and perhaps trial and error. The 
fact that it is the half-power point of the fundamental component of the 
approximation may be of some help. Also, note that as t — > °o , <£(t) — > 
Ax exp (— fl I t I ). We might choose 9. that I>(t) and *(r) approach 
zero at the same rate. However, matching autocorrelation functions by 
means of their tails is not necessarily a desirable approach. 
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8.2 System Analysis and Synthesis 

The noise variance ratio of a linear discrete system for which the 
input noise autocorrelation function has been approximated by a linear 
combination of wide-sense Markov autocorrelation functions may be 
obtained by substituting (59) in (5). We have 

M 2 -IEZ WiWfA, exp ( -Q V T \i-j\). (76) 

i=0 J=0 v=l 



Now let 



Then 



a, = exp (-O.T). (77) 



M 2 = E ATS ZV.^xt. 1 *-"]. (78) 

i-=i L »'=o i=o J 

The expression in brackets represents the noise variance ratio of the 
linear discrete system when the wth wide-sense Markov component of 
the noise is the input. Thus, it is clear that the principle of superposi- 
tion can be used to find the total noise variance ratio. Figs. 7, 8, 11, 12, 
and 13 may be applied to the wide-sense Markov components individu- 
ally. 

Use of the linear combination noise model will not be profitable in 
determining the optimum smoother of a class. There is a matrix inver- 
sion required [refer to (25)] which is more easily performed directly 
with the actual autocorrelation matrix. One should keep in mind that 
the noise variance ratio of a digital smoother is relatively insensitive to 
departures of the weighting function from the optimum. Hence a 
smoother optimized for the simple wide-sense Markov model may be 
satisfactory. 

The synthesis of a polynomial smoother based on the linear combina- 
tion noise model follows the method of Section VI, except that calcula- 
tion of 6 is somewhat more difficult, since oc must be calculated using 
(78) and the relevant plots. It should be noted that the estimate of e T 
may not be sufficiently accurate to justify the use of the linear combina- 
tion noise model. One should consider whether or not the simple wide- 
sense Markov model might be satisfactory. 

8.3 Noise Generation for Simulation 

Discrete stationary random noise of arbitrary autocorrelation func- 
tion <p(t) and variance a 2 may be approximately generated as a linear 
combination of independent, wide-sense Markov components. Let 
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Zi represent the ith sample of the approximating linear combination 

Zi = 2 b r Y ri , (79) 

where Y V i is the ith sample of the vih wide-sense Markov component. 
This i»th component is generated (refer to Section VII) as 

Y rl = X vl , (80) 



Y vi = a„F ril -i + Vl - ajX vi , (i > 1). (81) 

Care must be taken that the normalized uncorrelated random numbers 
A',., are generated in z similar but mutually uncorrelated sequences 
[X u \ t {Xu}, ..., [X,<] to ensure that the sequences \Y U }, {Y 2i }, 
j F„) are mutually uncorrelated. Note that each Y vi will liave zero mean 
and unit variance. 

To evaluate the coefficients b v , approximate the autocorrelation 
function of the arbitrary random noise process by a linear combination 
of wide-sense Markov components. Thus, from (59) and (77) we obtain 

*( \i - 3 I T) ^k\i-j\T) = JT A,**,}'-". (82) 

Now since the Z process is stationary 

z 

*(\i-j\T) = COv(Zi > Zi) 



var (Zi) °- 

tbiaj^ 



(83) 



We have set var (Zi) = a 2 since the arbitrary process and its approxi- 
mation must be matched in variance. Now, equating terms of (82) and 
(83), we obtain 

b v = aVT v . (84) 

Thus, 

Z, = aJ2VT v Y vi . (85) 

IX. GLOSSARY OF SYMBOLS 

A v = Coefficient in approximation of power density spectrum or 
autocorrelation function by linear combination of wide-sense 
Markov components 
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B = tlT s = noise-smoother bandwidth ratio for fundamental wide- 
sense Markov component 
C = output signal of smoother 

E = expected value operator = I dF 

v— 00 

f = dF = probability density function 

M = matrix of moments of weighting coefficients 

M q = 2 i"T''Wi = (jth moment of weighting function 

1=0 

n = present sample 

N = number of samples operated on by nonrecursive discrete 

smoother 
p = order of smoother 
P = autocorrelation matrix 
r = degree of smoother 
It = total input signal to smoother 
R m -i = total input signal evaluated at t = (m — i)T 
R = desired component of input signal to smoother 
]} (p) = pth derivative of desired component of input signal 
R = polynomial approximation to desired component of input signal 
R = noise component of input signal 
S(u) = power density spectrum (normalized sense, 

— / S(u) da = 1 ) ; Fourier transform of <p(t). 

i,r = time variables (seconds) 

T = sampling interval (seconds) 

T s = smoothing interval (seconds) 

Wi = weighting function of digital filter evaluated at t = iT 

X(t) = white noise process 

Y(t) = wide-sense Markov noise process 

Z(l) = general noise process 

Z(t) = approximation to general noise process 

a = exp( — 07') = intersample correlation for fundamental wide- 
sense Markov component 

a v = exp(— Q.T) = intersample correlation for t-th wide-sense 
Markov component 

r = prediction time 
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t = total output error of smoother 

e T = truncation error 

ix = ad or = output-input standard deviation ratio 

a = noise variance 

trc" = output noise variance 

a R ~ = input noise variance 

<p(t) = autocorrelation function 

co = angular frequency variable (radians/sec) 

ft = bandwidth of fundamental wide-sense Markov power density 

spectrum (radians/sec) 
!2„ = bandwidth of ?'th wide-sense Markov power density spectrum 

component (radians/sec) 
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APPENDIX 

Asymptotic Behavior of Smoothers 

We will consider the behavior of n as N —* °° with fi and T s fixed. 
We shall first find the limits of two expressions which will be needed in 
finding the limits of the larger noise variance ratio expressions: 

lim a aN+b = limexp[-flT(cW + &)] = lim exp \ - UT s (aN + b) l 

(86) 
.. f B(aN + b)l ( m 

= hmexp - — — = exp (-aB), 

JV-»oo |_ 1\ — 1 J 

1 -exp[-B/(iV- 1)] 



lim JV(1 — a) = lim 



-lim^^'-^-^-B. 

N^oo (N — l) 2 



l/N 

(87) 
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A.l Optimum Smoother — 0th Order, 0th Degree 

We have, using (33) 

i- 2 v 1 + a v 1 + a 2 , QQ v 

lim n = lim — j—z — — = lim -=-r- r — — — = . (.»»; 

A--.* K -<» N - (N - 2) a A--» N(l - a) + 2a B 4- 2 

A.2 Optimum Smoother — isi Order, 1st Degree 

We have, using (35) and (36) : 

.. T 2 2 _ ,. 12(l+q)[iV(l -a) - 1 + «] 

SS S M " i~ [iV(l - a) + 2«][JV(1 - «) + 1 + 3a] + 4a 

24J3 _ 24£ 

" (B 4- 2) (5 4- 4) 4- 4 ~ £ 2 + 65 + 12 ' 

A.3 Zeroth Order Cascaded Simple Averages Smoother 

Using (44) we have 

r 2 ,. / l + « , 2a(a Af - 1) 

lim u = lim S^tt^ r -f- t-KT/-, \v> 

at-oo^ *-*« (iV(l - a) [JV(1 - a)] 2 , 



= 1 + -S t ex P (-*) - 1] = -5 to C-B) +5-1]. 



A.4 First Order Cascaded Simple Averages Smoother 

We have, using (45), 
limT fl V=lim2(4.5) 2 (^- 1 Y 

J 14-q q(a W - 2q 2 » /3 - / /3 + 2) 

\3A/(l-a) ' [iV(l-a)] 2 (j[ 

= ^ |-exp (-£) + 2exp (-2B/3) 

4-exp(-£/3) +|JB-2J. 

A.5 Second Order Cascaded Simple Averages Smoother 
Using (46) we have 
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S^-S^f-rJfcfe 



+ 





a(a N - 2a w " - a™' 3 + 2a m + 3a™ - 3) 



[*(!-«)]' (92) 

= ^ jexp (-5) - 2 exp (-55/6) - exp (-25/3) 

+ 2 exp ( -5/2) + 3 exp -5/3) + 1 5 - 3} . 

REFERENCES 

1. Blackman, R. B., Smoothing and Prediction of Time Series by Cascaded 

Simple Averages, 1960 IRE Convention Record, 8, Part 2, March 21-24, 
pp. 47-54. Also published as Bell System Monograph 3678. 

2. Lees, A. B., Interpolation and Extrapolation of Sampled Data, IRE Trans. 

on Information Theory, IT-2, 1, March, 1956, pp. 12-17. 

3. Johnson, K. R., Optimum, Linear, Discrete Filtering of Signals Containing a 

Nonrandom Component, IRE Trans, on Information Theory, IT-2, 2, 
June, 1956, pp. 49-55. 

4. Blum, M., An Extension of the Minimum Mean Square Prediction Theory for 

Sampled Input Signals, IRE Trans, on Information Theory, IT-2, 3, 
September, 1956, pp. 176-184. 

5. Blum, M., On the Mean Square Noise Power of an Optimum Linear Discrete 

Filter Operating on Polynomial Plus White Noise Input, IRE Trans, on 
Information Theory, IT-3, 4, December, 1957, pp. 225-231. 

6. Doob, J. L., Stochastic Processes, Wiley, New York, 1953, p. 90. 

7. Doob, J. L., op. cit., pp. 233-234. 

8. Blackman, R. B., Linear Data Smoothing and Prediction in Theory and Prac- 

tice, to be published. 

9. Ragazzini, J. R., and Franklin, G. F., Sampled-dala Control Systems, McGraw- 

Hill, New York, 1958. 

10. Hamming, R. W., Numerical Methods for Scientists and Engineers, McGraw- 

Hill, New York, 1962, pp. 143-152. 

11. Blackman, R. B., Linear Data Smoothing and Prediction in Theory and Prac- 

tice. 

12. Taussky, Olga, and Todd, J., Generation and Testing of Pseudo-Random 

Numbers, Symposium on Monte Carlo Methods, Wiley, New York, 1956, pp. 
15-28. 

13. Chartres, B. A., An Exact Method of Generating Random Normal Deviates, 

T. R. No. 5, Contract No. DA-30-069-SC-78130, Division of Applied Mathe- 
matics, Brown University, March, 1959. 

14. Laning, J. H., and Battin, R. II., Random Processes in Automatic Control, 

McGraw-Hill, New York, 1956, pp. 381-394. 

15. Lee, Y. W., Statistical Theory of Communication, Wiley, New York, 1960, pp. 

460-469. 



